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Abstract We restudy the phase diagram of the 2D-Ising model with competing interactions Ji on nearest 
neighbour and J2 on next-nearest neighbour bonds via Monte-Carlo simulations. We present the finite 
temperature phase diagram and introduce computational methods which allow us to calculate transition 
temperatures close to the critical point at J2 = Ji/2. Further on we investigate the character of the 
different phase boundaries and find that the transition is weakly first order for moderate J2 > Ji/2. 

PACS. 05.50.+q Lattice theory and statistics including Ising, Potts models, etc - 75.10.Hk Classical spin 
models - 64.60.De Statistical mechanics of model systems (Ising model, Potts model, field-theory models, 
Monte Carlo techniques, etc) 



1 Introduction 

The search for exotic groundstates in two-dimensional frus- 
trated quantum spin systems is a topic of intense research 
(see, e.g., Refs. [1,2] for recent reviews). In this context, 
the antiferromagnetic J1-J2 spin-1/2 Heisenberg model on 
the square lattice has been intensively studied during the 
past two decades. Nevertheless, the nature of an interme- 
diate non-magnetic phase around J2 « Ji/2 has remained 
under debate until recently (see Refs. [1,3] and references 
therein). Among the possible approaches to the antiferro- 
magnetic J1-J2 Heisenberg model, quantum Monte-Carlo 
(QMC) simulations suffer from a severe sign problem in 
the region J 2 ~ Ji/2. Other approaches include perturba- 
tion theory around the Ising limit [4,5]. A closely related 
model is given by hard-core bosons on the square lattice 
with nearest neighbour (NN) and next-nearest neighbour 
(NNN) hopping and repulsion terms [6-9]. In this case, 
there is no sign problem such that QMC simulations are 
possible in principle [6-9]. However, simple QMC algo- 
rithms suffer freezing problems in the intermediate regime 
at J2 = Ji/2 which can be traced to a groundstate degen- 
eracy of the Ising limit. This motivated us to perform a 
model study by solving the related freezing problems in 
Monte-Carlo (MC) simulations of the Ising model. 

Investigations of the two-dimensional J\-J% Ising model 
on the square lattice have an even longer history than of 
the corresponding Heisenberg model, including in partic- 
ular MC simulations (see, e.g., [10-16]). Nevertheless, cer- 
tain issues have remained controversial also in the case 
of the Ising model, in particular the nature of the finite- 
temperature phase transition for J2 > Ji/2: MC simula- 
tions [10, 13, 15] and an investigation of the Fisher zeros 
of the partition function [17] have suggested non-universal 



critical exponents, whereas a variational approach [18] and 
a differential operator technique [19] predict a first-order 
transition for Ji/2 < J% < J\. In this paper we resolve 
this issue in favour of a weak first-order transition at least 
for not too large J2 > Ji/2 by providing substantially 
improved MC results for the phase diagram. 

This paper is organized as follows: in section 11.11 we 
introduce the J\-J% Ising model on the square lattice and 
discuss its T = groundstates. The MC simulation meth- 
ods are described in section [2] and results are presented 
in section [3J We conclude with a summary and outlook in 
section |4j 



1.1 Model 

We study the classical Ising model with competing anti- 
ferromagnetic interactions J\ on the NN bonds and J2 on 
the NNN bonds (Ji > 0): 
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We will study square lattices of linear extent L with peri- 
odic boundary conditions. 

For J2 = the groundstate of (p} is the known anti- 
ferromagnetic solution which is a Neel ordered lattice (see 
Fig. CD left) with energy E = -2J X N (N = L 2 is the 
number of sites). Switching on the repulsive interaction 
on the diagonal (NNN) bonds of the square lattice yields 
an increase of the groundstate energy for the Neel-ordered 
state: 
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Figure 1. Sketched are both ordered phases and a third 
groundstate configuration at J2 = Ji/2 (middle). The total 
degeneracy of the groundstate at J2 = Ji/2 is of order 2 L+1 . 
Shaded areas mark flipped lines. 



For large J2 the system orders in the collinear (or su- 
perantiferromagnetic) phase (see Fig. [TJ right) where all 
diagonal bonds are antiferromagnetic, while half of the 
NN bonds is antiferromagnetic and the other half is ferro- 
magnetic. Thus, the energy in this state depends only on 
J 2 - 



^Coll 



= -2NJ 2 



(3) 



The critical point separating these two phases lies at J2 = 
Ji/2, where the transition temperature is suppressed to 
T = 0. At this point the groundstate is highly degen- 
erate. More precisely, there are 2 L+1 — 2 groundstates 
which can be obtained as follows: flipping a line of an- 
tiparallel spins costs no energy at J 2 = Ji/2. Thus, one 
can generate almost 2 i+1 groundstates from a Neel state 
by flipping either the L horizontal or the L vertical lines 
independently (the second Neel state can be reached in 
both ways). In particular, one can reach a collinear state 
through a series of intermediate groundstates by L/2 line 
flips, as sketched in Fig. [TJ Accordingly, close to the criti- 
cal point J2 = Ji/2 the energy landscape is characterized 
by many local minima which are separated by large energy 
barriers. Therefore, MC simulations using only single-spin 
flips have problems to reach the configuration with global 
minimum energy. To solve this problem and to obtain 
transition temperatures also in the vicinity of J2 = Ji/2 
we implemented improved MC algorithms which we will 
discuss in the next section. 



2 Methods 

2.1 Computational methods 

We first tried to simulate the model with competing in- 
teractions using conventional single-spin flip MC simula- 
tions [16]. However, near the critical point this algorithm 
does not provide proper results and suffers severe freezing 
problems in the proximity of the critical temperature and 
below. To overcome these problems we implemented a par- 
allel tempering algorithm [20-23] : a number of simulations 



with the same set of parameters (Ji, J2, L) but varying 
temperatures are simulated simultaneously. After a suffi- 
ciently large number of sweeps over the complete lattice 
an additional MC step proposes a configuration swap be- 
tween neighbouring simulations with an acceptance rate 
p(i,i + l): 
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We have chosen the temperatures Ti logarithmically with 
a density maximum near the estimated transition temper- 
ature. A selfadjusting temperature set [23] was not neces- 
sary for our purposes. To manage the additional computa- 
tional workload of the parallel tempering algorithm we im- 
plemented a parallelisation of the MC code via OpenMP 
[24,25]. To assure the independence of the simulations we 
used the SPRNG (version 4.0) [26] lagged Fibonacci gener- 
ator for producing independent random number streams 
for each simulation. We have checked in some samples that 
our MC results do not change if we replace the random 
number generator by the Mersenne Twister algorithm [27] , 
as implemented in the BOOST libraries version 1.33.1. 

Another way to improve the MC simulations is to al- 
low not only single-spin updates but also flipping whole 
lines of spins in an additional MC step (as sketched in 
Fig. [T|) . This method simulates directly the transition be- 
tween degenerate groundstates at J2 = Ji/2 and helps to 
minimize statistical errors for lattice sizes up to 50 x 50. 
For larger lattices the probability to flip a whole line de- 
creases rapidly for ratios J 2 7^ ^i/2 while the simulation 
time to calculate them increases with L. We also tried 
cluster updates [28] but in the case of competing interac- 
tions this method does not help: when one approaches the 
critical temperature, the clusters extend over the whole 
lattice and therefore do not support the ordering process. 

The parallel tempering algorithm gives rise to correla- 
tions between different temperature points in a simulation 
which can affect the independency of the calculated data. 
Therefore, the meanvalues and errorbars of the shown ob- 
servables are derived from at least 10 independent MC 
runs. 



2.2 Order parameters 

To distinguish the two ordered and the disordered phases 
we use the respective structure factors 



(Xi —Xj ) 



(SiSj) 



(5) 



where q is a vector in momentum space and indicates the 
different magnetic phases: 

q = (0, 0) — > ferromagnetic order, 
q = (0, 7r), (it, 0) — > collinear order, (6) 
q = (%,tt) — > Neel order. 
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Figure 2. Binder cumulants depending on temperature for 
Ji — O.6J1 and different lattice sizes. In the inset three cumu- 
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Figure 3. Plotted are Tc's over the ratio J2/J1. Data were 



, , , , , , ., , , m , . , , . . m produced by parallel temperinq MC simulations. Near the crit- 

lants are plotted with errorbars. 1 he intersection area gives 1 c . , . , T m 1 i- ■ 17 r> 1 11 1 

ical point J2 = Ji/2 additional line flip updates led to a better 
with small error. c ' 

convergence. 



We identify M(q) — \J S(q)/N as our order parameter. 
The Neel order parameter can be calculated as staggered 
magnetization due to a simple sublattice rotation. Using 
the order of the collinear phase we can also calculate the 
associated order parameter via the column or line index of 
the underlying lattice. To find the transition temperature 
we calculate the fourth order Binder cumulant [16,29,30] 
for different system sizes L: 
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(7) 



For large enough L they will meet in a single point at 
Tc (for an example see Fig. [2]). To analyze the charac- 
ter of the phase transition we calculated the specific heat 
and recorded time series of the energies to set up his- 
tograms [31-33]. Furthermore we calculated the temper- 
ature derivative of the Binder cumulant which is related 
with the critical exponent v [16] to study the critical be- 
haviour of the system: 



dT 



(8) 
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To have a closer look on the critical point J2 = Ji/2 we 
analyzed the peaks of the specific heat via polynomial fit- 
ting. 



3 Results 

We calculated the critical temperatures for various ratios 
of J2/J1 especially close to the critical point (Fig. [3|). Our 
data is in good agreement with MC results from Lan- 
dau and Binder for J 2 > 0.6 Ji [13] and J 2 < 0.4 Ji [11]. 
In addition, our improved MC algorithm with the above 
described parallel tempering mode and line flip updates 
enabled us to obtain results much closer to J2 = Ji/2 



(| J2/Ji-0.5| = 0.005) for system sizes up to N = 250 000. 

To investigate the character of the phase transition we 
calculated the specific heat and energy histograms for dif- 
ferent lattice sizes. Looking at the peaks in the specific 
heat for different J2/J1 suggests different character of the 
phase transition. For J2 < Ji/2 (Fig. [Ha)) we have a 
slowly emerging peak with growing system size. Indeed, 
for an Ising like transition we expect only a logarithmic 
divergence in the specific heat [34] and a critical expo- 
nent v = 1 [35]. To verify the Ising like behaviour we 
had a closer look at the derivative of the Binder cumu- 
lant L*4 near the critical temperature. We expect a lin- 
ear correspondence in L for equation ([8]). An analysis of 
our results at Ji — 0.3Ji yields v = 0.99(1) which is in 
good agreement with the expected value v — 1. On the 
other hand for J2 > Ji/2 and in particular for the case 
J2 = 0.6 Ji shown in Fig. HJb) there are different opinions 
about the character of the phase transition. Some authors 
assume a continuous phase transition with non universal 
exponents [10,13,15,17] whereas other authors find a first 
order phase transition [18,19]. Note that our results for 
the specific heat (Fig. IHL)) differ quantitatively from the 
approximate results of [18]. In particular, the maximum 
of the specific heat continues to diverge for growing L. 
We estimated the area under the peak in dependence on 
the lattice size and find it to converge towards a constant 
value for large enough L. This indicates a S-peak struc- 
ture for the specific heat in the thermodynamic limit as 
expected for a first order phase transition. We should nev- 
ertheless mention that the value of the specific heat does 
not follow the finite size scaling law expected for first order 
transitions [31] very well, which may be due to crossover 
phenomena. 

As further verification of our characterisation of the 
order of the phase transition we computed energy his- 
tograms. First, we extracted a system-size dependent Tc 
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Figure 4. Comparison of specific heats for J2 < Ji/2 (a) and 
J2 > Ji/2 (b). The peaks emerge at the critical temperature 
obtained from the fourth order cumulant calculations. Note the 
different magnitudes of the divergence. Errorbars are omitted 
if they are of the same order as the linewidth. 
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Figure 5. Energy histogram for J2 = 0.3Ji (binsize — 
10 _4 Ji) and J2 = O.6J1 (binsize — 6 ■ 10~ 3 Ji). For each lat- 
tice size L the histogram is plotted at the critical temperature 
Tc(L). The slowly emerging double peak structure in panel (b) 
in comparison to the clearly single peaked structure in panel 
(a) indicates a weak first order transition for J2 > Ji/2. 



from the maximum of the specific heat. Then we recorded 
the time series of the energy for each system size and the 
associated Tc- Fig. [5] shows the resulting energy distribu- 
tion for J 2 = 0.3Ji and J 2 = O.6J1. For J 2 = 0.3Ji we 
find a single pealQ, consistent with a conventional second 
order phase transition. By contrast, for J 2 = 0.6 Ji a dou- 
ble peak structure emerges for sufficiently large system 
size, see Fig. [5jb). Such a double peak structure is char- 
acteristic for a first order transition [31-33]. The fact that 
the double peak structure emerges only for large system 
sizes shows that the transition for J 2 = O.6J1 is a weak 
first order one. Since we have observed similar behavior 
for nearby ratios of J 2 /Ji, we believe the transition for 
J2 > Ji/2 to be of first order, at least for not too large 
values of J 2 . 

Fig. [6] shows the specific heat exactly at the critical 
point J 2 = Ji/2 for different system sizes. The purpose of 
this computation is to verify if a direct phase transition 
from Neel to collinear order is possible at finite tempera- 
ture, or if any other finite-temperature phase could exist 
in this region. Previous publications [11, 18] have shown 
curves for the specific heat at J 2 = Ji/2 with a rounded 
peak, but we are not aware of any systematic finite-size 
analysis. In our results, we observe that the peaks of the 
specific heat are moving to lower temperatures for increas- 
ing system sizes, suggesting that the transition tempera- 
ture is suppressed to Tc = in the thermodynamic limit. 
To further substantiate this conclusion, the inset of Fig. [5] 
shows the peak positions (Tc) with respect to the inverse 
lattice length \/L. For the given lattice sizes we find a 
power law behaviour for Tc and l/L. This is consistent 
with Tc — in the thermodynamic limit. 



1 The histogram shown in Fig. [5ja) is very close to a single 
Gaussian curve, as is known for second order phase transitions 
[32]. 



0.01 1/L 0.02 0.03 




T/Jj 

Figure 6. Specific heat at J2 = 0.5Ji for different lat- 
tice sizes with errorbars at every fourth data point. The 
peak moves towards lower temperature (here Tc/Ji = 
0.380(4), 0.345(4), 0.331(3), 0.320(2)) for growing system sizes. 
The inset shows Tc (obtained from a polynomial fitting of the 
maximum in the specific heat) versus the inverse lattice length 
in a log-log-scale (errorbars are smaller than the symbols) . The 
behaviour of Tc is consistent with a power law in L, indicating 
Tc = in the thermodynamic limit. 



4 Discussion 

In this paper, we combined single-spin flip, parallel tem- 
pering and line flip MC algorithms to enhance the effi- 
ciency of MC simulations of the frustrated square-lattice 
Ising model. These improvements enabled us to compute 
critical temperatures in the direct vicinity of the critical 
point J 2 = Ji/2 where the groundstate is highly degen- 
erate. For J 2 < Ji/2 there is a finite-temperature phase 
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transition into a Neel-ordered state. This transition be- 
longs to the two-dimensional Ising universality class. We 
believe that at the critical point J 2 = Ji/2 the transi- 
tion temperature is suppressed to zero. On the right hand 
side (J2 > Ji/2) there is again a finite-temperature phase 
transition into a phase with collincar order. 

At J 2 = O.6J1 we have observed a double peak struc- 
ture in the histogram of the energy. This identifies the 
phase transition as a weakly first order one. We believe 
this first order transition to be generic for moderate J2 > 
Ji/2. However, for larger J 2 it becomes increasingly diffi- 
cult to identify a double peak structure in the histogram 
of the energy. Therefore, the transition could in fact be- 
come a second order one for J 2 > J\. Note that a recent 
MC investigation [15] has concluded that the transition is 
second order for J 2 = Ji, although in [15] only smaller 
lattices have been considered than in the present work. 

In order to understand the possible nature of the tran- 
sition at large J2, it is instructive to consider the limit 
J\ = where one has two decoupled Ising models. The 
critical theory then consists of two conformal field theo- 
ries with central charge c = 1/2 each, i.e., a c = 1 theory. 
Thus, the point J\ = lies within the manifold of c = 1 
conformal field theories where continuously varying criti- 
cal exponents are possible in principle [36]. However, the 
coupling given by J\ turns out to be a relevant pertur- 
bation of the critical point such that this coupling should 
either lead to a first order transition or a fixed point with 
c < 1 [37] . According to the classification of minimal con- 
formal field theories with c < 1 [38] a possible second 
order phase transition at J 2 » J\ should have universal 
exponents. Therefore, a scenario with non-universal criti- 
cal behaviour [10,13,15,17] does not appear very plausible 
from a conformal field theory perspective either. Our MC 
results show that large crossover scales exist in the present 
model such that very big lattices would be needed for a 
reliable numerical determination of the nature of the tran- 
sition for J2 > J\. 

A similar phase diagram as in the Ising model is also 
found for the classical J1-J2 X-Y [39] and Heisenberg 
models [40]. However, we believe that the nature of the 
phase transition for J 2 > J\/2 is a different issue in these 
two models due to the different symmetries in spin space. 

After having solved the freezing problems in the Ising 
model, we are now about to introduce hopping terms and 
to study finite-temperature properties of hard-core bosons 
on the square lattice using parallel tempering QMC sim- 
ulations [41]. 

We acknowledge financial support by the Deutsche Forschungs- 
gemeinschaft under grant No. HO 2325/4-1 and through SFB602. 
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